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I This work considers the problem of locating a single source from noisy range measurements to 

(~| , a set of nodes in a wireless sensor network. We propose two new techniques that we designate as 

. Source LocaUzation with Nuclear Norm (SLNN) and Source Localization with ^i-norm (SL-£i), which 

extend to arbitrary real dimensions, including 3D, our prior work on 2D source localization formulated 

in the complex plane. Broadly, our approach is based on formulating a Maximum-Likelihood (ML) 

\ estimation problem for the source position, and then using convex relaxation techniques to obtain a 

lO i semidefinite program (SDP) that can be globally and efficiently solved. SLNN directly approximates 

l> ■ 

\^ • the Gaussian ML solution, and the relaxation is shown to be tighter than in other methods in the same 

class. We present an analysis of the convexity properties of the constraint set for the 2D complex 
version of SLNN (SLCP) to justify the observed tightness of the relaxation. In terms of global accuracy 
of localization, SLNN outperforms state-of-the-art optimization-based methods with either iterative or 
. ; closed-form formulations. We propose the SL-^i algorithm to address the Laplacian noise case, which 

' models the presence of outliers in range measurements. We overcome the nondifferentiability of the 



Laplacian likelihood function by rewriting the ML problem as an exact weighted version of the Gaussian 
case, and compare two solution strategies. One of them is iterative, based on block coordinate descent, and 
uses SLNN as a subprocessing block. The other, attaining only slightly worse performance, is noniterative 
and based on an SDP relaxation of the weighted ML problem. 
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I. Introduction 

Locating a source from range measurements to a set of known reference points (anchors) is a classic 
problem in many engineering applications (e.g., radar, sonar, GPS), and has received a great deal of 
attention over the years. Recently, source localization from range measurements has been intensively 
examined in the context of wireless sensor networks (WSN), where ranges estimated from times of arrival, 
or from surrogates such as received signal strength, are somewhat unreliable due to the complexity of 
many WSN propagation environments (e.g., indoor settings with few unobstructed line-of-sight paths). 

Spatial information per se, or as georeference to other sensor measurements, is crucial in WSN 
applications and warrants investigation into suitable localization algorithms. While many approaches to 
source localization based on classical triangulation or heuristic criteria can be found in the WSN literature 
m, m, our primary focus is on optimization-based methods formally derived from the likelihood function 
of observations, or related cost functions S-IH. By doing so, we expect to take advantage of the 
optimality properties of maximum likelihood (ML) estimates to improve the robustness to perturbations 
in range measurements. We do not consider alternative/complementary measurements such as angles of 
arrival or time differences of arrival. We also assume cooperative localization scenarios where absolute 
ranges, as opposed to range differences to a reference sensor, are measured. These can be obtained either 
by synchronizing clocks and transmitting waveforms from the source at known times (beacon mode), or 
by initiating the transmission at a reference sensor and measuring the round trip time to the source and 
back (transponder mode). 

Centralized ML algorithms for range -based source localization, which require the transmission of the 
full data set to a fusion node for processing, are proposed in Q, HI under Gaussian noise and in 
|[9l under Laplacian noise. These resort to semidefinite relaxation (SDR) to alleviate the problem of 
algorithmic convergence to undesirable local maxima of the likelihood function. A related alternative 
approach proposed in |4| solves a constrained least-squares (LS) problem using squared range (SR) 
measurements, subject to a quadratic constraint. This was shown to outperform, on average, the ML SDR 
approach of lEl whose relaxed solutions sometimes fail to produce meaningful source position vectors 
(rank one solutions). Another approach, proposed in [SI, approximates the ML solution via second-order 
cone programming and a low-dimensional search. 

Distributed algorithms for wireless sensor nodes, where the source location is iteratively determined 
through in-network processing at individual nodes and communication between neighbours, are also being 
very actively pursued |fT0l - |[T2ll . These techniques, however, are not the focus of our work. We also note 
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that source localization can be viewed as a special instance of sensor network localization, where the 
positions of several sources/sensors are simultaneously determined from pairwise range measurements. 
Related algorithms based on semidefinite programming (SDP) have been developed for this class of 
problems |[T3i . and are relevant when there is significant uncertainty in anchor positions (see, e.g., 
llT4l for a similar SDP approach to source localization with anchor uncertainty using range differences). 

This paper develops an alternative to the source localization ML SDR method of |8]. We term this 
approach, originally proposed in fT] for 2D locaUzation under Gaussian noise. Source Localization in the 
Complex Plane (SLCP). Our relaxation for the nonconvex and nonsmooth likelihood function is tighter 
than the one presented in HI, in the sense that the relaxed solution will more often have (near) rank-1, 
as required to obtain target coordinates by factorization. SLCP also outperforms the SR-LS method of 
im, which iteratively solves a generalized trust-region subproblem and dispenses with factorization of 
rank-1 matrices, but undergoes some degradation with noisy measurements due to squaring of ranges in 
the cost function. The degradation of SR-LS becomes more severe in the presence of outliers ||9l, which 
commonly affect practical range measurement systems, e.g., when non-line-of-sight propagation occurs. 

This paper expands upon the results of |3 | in several ways: 

1) We extend the framework of SLCP from 2D localization, which relied on a formulation where 
target and anchor coordinates were represented as complex numbers, to arbitrary (real) dimensions. 
We term the new SDR method Source Localization with Nuclear Norm (SLNN), as this norm arises 
naturally in the cost function of our relaxed optimization problem. Similarly to SLCP, SLNN offers 
a tight relaxation in most problem instances, and retains a performance advantage over SR-LS. 

2) We provide a more complete analysis of the accuracy properties of SLCP, whose success in 
providing tight relaxations relies on certain parametrically defined sets in being nearly convex. 
We discuss the convexity of the sets and how to trace the convex hull for any of them, from which 
convexity can be empirically assessed. For three-anchor scenarios we also examine a search-based 
alternative to SVD decomposition to extract the source coordinates from the solution of SLCP (a 
positive semidefinite matrix with near rank-1). 

3) In 191 a modification of SLCP, termed SL-£i, was introduced for ML source localization under 
Laplacian noise. This makes the algorithm robust to outlier measurements, a property that was 
observed in simulation even for non-Laplacian range errors. In this paper we provide a conceptually 
similar extension for source locaUzation beyond 2D, consisting of a reformulation of the nondif- 
ferentiable log-likelihood function for Laplacian noise as a reweighted version of the Gaussian 
log-likelihood. We propose both single convex formulations and a simpler iterative optimization 
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algorithm, which repeatedly solves weighted SLNN problems followed by weight refinement. 

The organization of the paper is as follows. In SectionUll we formulate the ML source location problem 
under Gaussian or Laplacian noise. In Section |lll] we derive the SLCP algorithm for 2D localization, 
we analyze the geometry of the associated optimization problem and the tightness of the relaxation 
(Subsection IIII-AI) . and we propose criteria for factorizing the SDR solution to recover the source 
coordinates (Subsection IIII-BI ). In Section |IV] we derive the SLNN algorithm, which extends SLCP to 
3 (and higher) dimensions, and propose an iterative version of this algorithm that can handle Laplacian 
noise (Section lIV-Ab . Section |V] illustrates the performances of the algorithms in simulation. Finally, 
conclusions are drawn in Section IVll 

Throughout, both scalars and individual position vectors are represented by lowercase letters. Other 
vectors and matrices are denoted by boldface lowercase and uppercase letters, respectively. Individual 
components of matrix X are written as Xij and those of vector x as Xi (the same notation would be 
used for a hypothetical position vector Xj, but the distinction between both should be clear from context). 
The superscript T (H) denotes the transpose (hermitian) of the given real (complex) vector or matrix, 
(•, •) denotes the inner product of two vectors, and tr(-) denotes the trace of a matrix. For symmetric 
matrix X, X ^ means that X is positive semidefinite. We denote the Frobenius norm of matrix X as 
||X||ir = ^tr(X-^X) and its nuclear norm as ||X||Ar = tr((X^X)2). Below, is the m x m identity 
matrix and 1^ is the vector of m ones. The convex hull of set S is denoted by co(5). 

II. Problem formulation 

Let X G M" be the unknown source position, Oj € M", i = I, ..,m be known sensor positions (anchors), 
and Tj = — ajll be the measured range between the source and the i-th anchor, where Wi denotes a 
noise term with standard deviation a. Under i.i.d. Gaussian or Laplacian noise maximizing the likelihood 
of observations for the source localization problem is equivalent to 

minimize YlT^-i Ilk — aJI^ — rf|''. 

'-^ ' (1) 

X 

We will derive the SLCP/SLNN algorithms to (approximately) solve ([T]) under Gaussian noise (jj = 1, 
q = 2), whereas SL-£i will solve it under Laplacian noise (p = I, q = 1). The case (p = 2, q = 2) is 
also of interest and corresponds to the cost function used in the SR-LS algorithm of Pl, which is used 
to benchmark our algorithms. Note that the cost function for SR-LS is not a likelihood function, and it 
arises out of mathematical convenience. 
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Fig. 1: Geometrical interpretation of terms in the source localization cost function ([T]) for p = 1, q = 2. 

The main difficulties of solving ([T) lie in the fact that this cost function is, in general, nonconvex 
and multimodal. For q = 1 it is also nondifferentiable, which poses additional challenges. We address 
the nonconvexity and multimodality of the cost function in Sections Hn] and |IV] by developing convex 
relaxations that turn out to be tight in most problem instances, thus providing a very good approximation 
to the true source location. If necessary, the source coordinates can be further refined by iteratively 
minimizing ^ starting from the relaxed solution. Refer to |9] for one such iterative refinement approach 
based on the Majorization-Minimization algorithm. We address the nondifferentiability of ([T|l for q = 1 
in Section IIV-AI by rewriting it as a weighted version of the case q = 2, where the weights themselves 
become optimization variables. 

III. Source Localization in 2D: SLCP 

For p = 1, g = 2 we view each term in ([T]) as the squared distance between two circles centered on 
aj, one with radius ||x — aj||, and the other with radius rj (see Figure [T|). This term can be replaced 
by the squared norm of the difference between the position vector x and its closest point on the circle 
{y ^ M? : \\y — ai\\ = r.i} , which we denote by y,. Problem ([T]) can then be equivalently expressed as (a 
formal proof of equivalence is provided in |9l) 

minimize YlTLi Ik - 

X, Vi (2) 

subject to \\yi — ai\\ = ri i = 1, . . . ,m. 
If we fix yi, the solution of ^ with respect to x is an unconstrained optimization problem whose 
solution is readily obtained as the center of mass of the constellation x = ^ Sl^i Vi- Moreover, in 2D 
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the constraints of Q can be compactly described in the complex plane, yielding 

minimize H^^l^l^y-yp 

y,9 (3) 

subject to y = a + HO, 



where a 



ai 



T 



G C™ holds the anchor coordinates, expressed as complex numbers, R = 

r 1 ^ 

diag(ri, . . . ,rm) G M™^"", and = e^'f'^ . . . e-''^" £ C™. The complex representation makes it 



simple to impose unit magnitude constraints on the elements of 9, and later relax them to obtain an SDR. 
Expanding the objective function and deleting constant terms yields the quadratic constrained problem 

minimize 2Re(c^ 9) - ^9^rr'^9 

e (4) 

subject to = 1, 
where r = Rl^ and c = R(Im - ^lml^)a. 

To proceed we now wish to replace Re{c^9) in ^ with —\c^9\, which is readily written as a function 
of a quadratic form in 9 and then relaxed in the same way as the second term in the objective function. 
To this end, first note that if 9 is replaced with 9e^'^ neither the second term in the objective function 
of ^ nor the constraints change for any angle 7. By proper choice of 7 the complex number c^9 may 
be rotated to the (negative) real axis for any feasible 9, such that Re{c^ 9e^'^') = —\c^9\, thus reducing 
the value of the objective function relative to other values of 7. This implies that any optimal solution 
of ^ will satisfy Re(c^0) = —\c^9\, which justifies replacing Re(-) with — | • | in the cost function. It 
should be kept in mind, however, that once a solution 9 to the modified optimization problem is obtained 
it should be rotated to obtain the actual vector of phases 9e^'^ such that Re(c^0e-^^) = —\c^9\. 

Now the modified problem is equivalently written as 



maximize 2^ti{cc^99") + ^ti{rr^99^) 

9 (5) 
subject to = 1, 

and following standard manipulations we introduce the new variable $ = 99^ and an associated 
(nonconvex) constraint rank($) = 1. Finally, a SDR formulation of SLCP is obtained by introducing the 
hypograph variable t such that < t < 2Y/tr(cc^$) and dropping the rank constraint 

maximize t + — tr(rr^$) 

(6) 

subject to * ^ 0, (pa = 1, 4c^$c > t'^. 
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TABLE I: Summary of the SLCP algorithm 



1) Given the anchor positions and range measurements, solve the SDR ([SJl 

2) Compute a rank-1 approximation of the SDR solution as $ « 00^ 

3) Compute a rotation angle 7 such that Re(c^0e-'^) = — |c^0| in Q 

4) Obtain the vector of circle projections y = a + Rffe-''' 

5) Estimate the source position as the centroid x = ^l^y 

Remark that the solution of ([6]) is a positive semidefinite matrix, which should have a clearly dominant 
eigenvalue in problem instances where the SDR is an accurate approximation to the initial problem In 
such cases $ ^ Aiuiuf^, where Ai is the highest eigenvalue of $ and ui the corresponding eigenvector, 
and the vector of complex phases is estimated as = y/Xl\ii |15 |. An alternative approach for computing 
6 is examined in Section IIII-BI Table U summarizes the SLCP algorithm. 

A. Tightness and Geometry of the Constraint Set in SLCP 

The source localization problem prior to relaxation ^ can be written as 

maximize 2-v/u + —v 

V m 



The objective function in ^ is concave with respect to u and v, and the optimization problem would be 
convex if the set S, over which this function should be maximized, were convex. Then, the SDR used 
in SLCP Q would always find a rank-1 solution from which the vector of phases 6 would readily 
follow by factorization. In practice it was found that, even for a moderate number of anchors, the set S 
is likely to have the required shape along part of its border, as discussed below, so that the SDR solution 
has indeed rank-1. We now examine some of the properties of S and the optimal solution. 

Given the separable form of the cost function ([7]) it is clear that, for fixed v, it can be maximized by 
choosing u as large as possible within S, and vice-versa. This implies the following property for the 
optimal points of ([7]): 

Property 1. The optimal points of ([7]) lie on the "upper right" boundary of set S, i.e., optimal points 
of ^ are maximal elements of S with respect to the standard cone Ml /ii6|/ . 



u, V 



(7) 



subject to {u, v) e S 



where 




(8) 
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Regarding the convexity properties of S, recall that the cost function of ^ was designed to be invariant 
to rotations of 6 so that, without loss of generality, the first element may be taken as unity. For m = 2 
anchors and 0i = 1, O2 = e^'^ we then have 

u = \c*i + C2e^''^p = |ci|^ + |c2|^ + 2|ci||c2| cos((/> + a) (9) 
V = \ri+ r2e>'^\^ = rl + rl + 2rir2 cos (10) 

where a = Zci — Zc2. Set S is an ellipse centered on (|cip + |c2p, rl + r"^), therefore clearly nonconvex. 
Given the definitions of c and r in for m > 2 anchors it is always possible to zero out elements 
3, . . . , m in these vectors if rs = . . . = r^, = in the diagonal of R, thus reverting to the case m = 2. 
In summary: 

Property 2. Depending on the specific range measurements, set S may be nonconvex for any number of 
anchors. 

In spite of the lack of convexity guarantees for S, our simulation results suggest that for m > 3 
anchors and typical range measurements this set usually does have a convex-like shape. Even when S 
is not convex all that is required for our SDR to provide a rank-1 solution is "local convexity" along 
the "upper right" boundary of S where the optimal point of ^ is known to be located. More formally, 
we require that the intersection of S with any supporting hyperplane defined by a normal direction with 
nonnegative components be a compact subset (a single point or a line segment) |fT6l . Figure |2] depicts 
some examples of S for different numbers of anchors and randomly generated c, r. Our practical test for 
local (non)convexity of S consists of tracing multiple supporting hyperplanes with nonnegative normal 
elements, and assessing whether any of them intersect S at two well-separated points. We build supporting 
hyperplanes not on S directly, which is a hard problem, but on the related relaxed convex set 

r = {(tr(cc^*), tr(rr^*)) : * E c™''", * ^ 0, <l)ii = l] . (11) 

Specifically, for a supporting hyperplane with normal (cos/3, sin/3), < /3 < ^, we determine an 
intersection point with T by solving the convex optimization problem 

maximize ((cos/3, sin/3), (tr(cc^$), tr(rr-^$))) 

* (12) 

subject to $ ^ 0, </>jj = 1, 
and setting the intersection point as n = tr(cc^$), v = tr(rr^$). This procedure is justified by the 
following result, proved in Appendix lAl 
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Fig. 2: The constraint set S for randomly generated 6 satisfying \9i\ = 1 and different numbers of 
randomly placed anchors. For each set the hypothesized convex hull, computed by relaxation of S, is 
also depicted. 



Lemma 1. For m < 3 anchors the sets S and T have the same set of supporting hyperplanes with 
nonnegative normal elements. Equivalently, in the relevant portion of its boundary T coincides with the 
convex hull of S. 

Although we only prove this result up to m = 3, the empirical evidence suggests that it is also valid for 
higher m, at least up to some maximum order (see Figure ll]). We leave this as a conjecture and apply the 
procedure for m > 3 as well, noting, however, that the case m = 3 has major practical significance as the 
minimum number of anchors that are necessary to recover a general 2D source position based on range 
measurements. We also conjecture that T is actually the convex hull of S, so (fT2l ) may be used to trace 
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the full boundary of co(5), and not just the portion where the supporting hyperplanes have nonnegative 
normal elements. This assumption is not required for our analysis, but was used for generating the set 
boundaries shown in Figure |2] 



B. Factorization of the SDR Solution 

The solution of the relaxed SLCP optimization problem Q is a positive semidefinite matrix, from 
which the vector of complex exponentials 6 is calculated by rank-1 factorization. The latter is needed 
to form the vector of circle projections y = a + HO (see ([3]l) and, ultimately, the source position vector 
as the centroid x = ^l^^y. The rank-1 factorization method advocated at the end of Section Hill is 
truncation of the eigenvalue decomposition of $ at the highest eigenvalue. Here we examine a more 
exact search-based alternative for the practically relevant case of m = 3 anchors, which will also be 
useful to assess the accuracy of the factorization based on eigenvalue truncation. 

For a given positive semidefinite matrix $ G Qmxm wish to find vector 6 G C™ satisfying 

minimize W^ — OO^Wj^ 

e (13) 

subject to = 1. 
The objective function in ([T3]) is expanded as 

||# - ee^fp = tr((* - ee^)^{^ - ee")) = ||*||| + \\ef -\i{i^"ee") - ti{ee^^). (U) 

Ignoring constant terms the optimization problem is equivalently reformulated as 

maximize O^^O 

e (15) 

subject to = 1. 

The cost function of dTSl ) is insensitive to a global rotation of all elements of by a common factor, 

r -iT 

hence for m = 3 anchors 6 can be written as = 1 e-'" and (fTSl) becomes 

maximize Re(^i2e^" + (/>23e^'^ + 0i3eJ("+^)). 
a, 5 

For fixed a the maximum is attained for 5 = —Z{(p23 + (pise^'^), yielding for ([T6l ) 

maximize Re((/>i2e-'") + |023 + ^ise-'"!- 
a 



(16) 



(17) 
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The solution to (fTTl ) is found by searching for the maximum value over the interval [0, 2tt). 

Referring to the definitions of the 2D sets 5 in ([8]l and T in (fTTl) . similar criteria to the above were 
considered for finding 9 such that the induced point in S is closest in Euclidean norm to the one induced 
by $ in T. However, the many-to-one nature of the mapping of onto points in S makes this formulation 
intrinsically ambiguous. 

IV. Source Localization in Higher Dimensions: SLNN 

To extend the approach used in SLCP to n > 2 dimensions, we write the circle/sphere equations in 
^ using an equivalent parametric form with real coordinates 

minimize \\x — ViW'^ 

x,yi,Ui (18) 



subject to yi = ai + riUi, 



h 



where x,yi,ai and Ui are now vectors in M", rather than complex scalars used in SLCP. In (fTSl ) Ui G 
is a unit-norm vector that plays the same role as the complex phase shift e^'^^ in SLCP. Equivalently, 



minimize 

XiVii'^i 



I Im,^ 





T 




T 




u{ 


subject to 


T 
Vm 




T 


+R 


T 



(19) 



Y A U 

where R = diag(ri, . . . , r^) as in For fixed yi, Ui ([T9l ) describes n uncoupled least-squares problems 
whose variables are the components of the source location vector x. The optimal solutions may be jointly 
written compactly as 



x^' = (1^1^)-4^Y = -l^Y. (20) 

m 

Replacing this back in ( fT9l ) to eliminate variable x the objective function becomes ||nY|||. = tr(Y^nY), 
where 11 = 1^ — ^Imlm ^ projection matrix (hence idempotent). Similarly to ©-(lU) we can now 
eliminate variable Y and the first set of equality constraints, expanding its definition in the objective 
function and ignoring constant terms to obtain 

minimize 2tr(C'^U) - :^tr(U'^rr'^U) 

U (21) 
subject to = 1, 

where C = RIIA and, as in (01), r = Rim- 



1 
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a ) Nuclear Norm Approximation: As in the complex formulation we wish to rewrite the first term in 
the objective function of ((2TI) in a form that is more amenable to SDR. In the optimization problem we thus 
replace U with the product UV, where V is an n x n orthogonal matrix such that V^V = VV^ = I„, 
yielding 



Note that, due to the orthogonality of V, each line of UV still has unit norm, so for any feasible U in (|2T]) 
UV is also feasible. Regarding (l22l ). V may be interpreted as an inner optimization variable that, for each 
candidate U, minimizes the value of the objective function. Noting that the second term in the objective 
function ^ does not depend on V, as tr(V^U^rr^UV) = tr(rr'^UVV'^U'^) = tr(rr^UU^), the 
inner optimization problem simply becomes 



This involves the minimization of a linear function on the set of orthogonal matrices, which resembles 
the known problem of minimizing a linear function of a vector v, say, (v, a), on the unit sphere ||v|p = 
v^v = 1. Invoking the KKT conditions |fT6l the latter problem is readily seen to yield the optimal cost 
— ||a||, attained at the point on the sphere along vector —a. One would therefore expect the solution 
of ( [23] ) to be — ||C^U||, involving some matrix norm of C^U. In Appendix IB] it is shown that this is 
indeed the case, and that the appropriate norm to consider is the nuclear norm, defined for matrix X as 
IIXIIat = tr((X^X)2), and equaling the sum of its singular values lITTl . The optimization problem (l22l ) 



minimize 2tr(C^UV) - itr(V^U^rr^UV) 
U,V 

subject to = 1, V-^V = I„. 



(22) 



minimize tr(C^UV) = (V, U^C) 



V 



(23) 



subject to V'^V = L 



is therefore equivalently rewritten as 



minimize 



2||C^U|U - ^tr(rr' 




U 



(24) 



subject to 



or 



maximize 2tr((C'^UU'^C) = ) + ^tr(rr'^UU^) 



U 



(25) 



subject to ||uj|| = 1. 
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TABLE II: Summary of the SLNN algorithm 

1) Given the anchor positions and range measurements, solve the SDR ( 127b 

2) Compute a rank-n approximation of the SDR solution as W ~ UU"^ 

3) Solve the inner optimization problem ([23} to get the rotation matrix V 

4) Obtain the matrix of sphere projections as Y = A + RUV 

5) Estimate the source position as the centroid of the rows of Y, a; = ^Y"^!™ 



We now introduce the variable W = UU^ and ignore the associated nonconvex constraint rank(W) = n 
to obtain the SDR 

maximize 2 tr((C^WC)^) + ^tr(rr'^W) 

W (26) 
subject to W ^ 0, Wii = 1. 
The objective function of (l26l ) is the sum of a concavq^ function of W with a linear term, and is therefore 



concave. The constraint set of ( |26l ) is convex, thus estabUshing that this is indeed a convex optimization 
problem. We express it in standard SDP form as 

maximize 2 tr(Z) + ^tr(rr-^W) 
W,Z 



subject to W ^ 0, Wi 



(27) 

C^WC z 

■ ^ 0, zy 0. 



Z In 



The equivalence between ( [26l ) and ( [27l ) is proved in Appendix |Bj 

Similarly to the complex 2D formulation, the solution of our SDR is a m x m matrix W that should 
have approximately rank n when the relaxation is tight. The matrix U of unit-norm vectors is obtained 
by SVD factorization of W |[T5l and, after accounting for the inner rotation of U, it is used to build the 
Hi and, ultimately, the source position vector x. Table |ll] summarizes the SLNN algorithm. 

A. Localization under Laplacian Noise: SL-ii 

When disturbances are Laplacian and i.i.d., thus heavier tailed than Gaussian, maximizing the likelihood 
amounts to solving ([T|) for p = g = 1, 

minimize i \ \\x — a,- 11 — rA. 

X 

'The first term is the composition of the linear map X — C^WC with tr(X2), which is known to be concave in X |16| . 
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The presence of |-| in each summation term of (l28l) . rather than (•)^, de-emphasizes the contributions of 
measurements rj corrupted by large noise values. The optimal point of (1281 ) is thus less biased by these 
outlier measurements than the cost function ([T|) for the Gaussian case p = I, q = 2. However, a major 
difficulty in solving (1281) is the fact that the cost function is not differentiable, making it less amenable to 
the types of analytic manipulations that we use to develop SDR. The strategy that we adopt to circumvent 
this difficulty parallels the one used in |9] for 2D sources, and as a key ingredient involves squaring the 
cost function of (1281 ) (which does not affect the location of extremal points), and then rewriting it as 

1=1 — — 

x,X (29) 
subject to Aj > 0, l^-^ = 1- 
The cost function is thus reduced to a weighted version of the more tractable Gaussian log-likelihood, 
where the real weighting coefficients Aj become optimization variables themselves. See |9| for a proof of 
this result (also |[T8l ). Now, the manipulations used earlier in Section |IV] for the development of SLNN 
can be replicated here to reformulate the problem as 

i=i — 

x,yi,Ui,X (30) 
subject to yi = ai + riUi, \\ui\\ = 1, Aj > 0, l^A = 1. 
For given yi, Ui, and A, (l30l ) has a least-squares cost function whose unconstrained optimal solution with 
respect to x is readily found in closed form from the first-order stationarity condition 

1=1 



(31) 



Substituting the optimal x in (l30l ). and using matrix notation, the cost function becomes tr(Y^HY), 
where H is the modified projector 





Ai 



1 



Em J_ 



J_ 

Ai 



L Ai ••• A„ 



(32) 



with A = diag(Ai, . . . , A^). 

1) Alternating directions (SL-£i AD): One possibility for iteratively solving (l30l ) is to use block 
coordinate descent, alternating between minimizing the expression with respect to {x, yi, Ui} for fixed 



November 30, 2011 



DRAFT 



14 



A and vice-versa. For fixed A tlie problem is 

minimize tr((A + RU)'^H(A + RU)) 

U (33) 
subject to ||ui|| = 1, 

which differs from the SLNN formulation only in the projector matrix S. It can therefore be similarly 
manipulated into a relaxed form that parallels (|27] ) 

maximize 2 tr(Z) + itr(rr^W) 
W,Z 

C^WC z 



subject to W ^ 0, Wii = I, 



(34) 



^0, Z ^ 0, 



with C = RHA, r = R 



_i_ 

Ai 



1 T 



, and K = ^. For the converse block coordinate descent 



step with fixed {x, yi, Ui} the problem is 

minimize -j^ 
A 

subject to Aj > 0, l^-^ = 1) 



(35) 



where Ki = |||x — aj|| — rj| = ||x — yj|| are constant in this subproblem. The solution, readily obtained 
from the first-order KKT conditions, is given by ||9l 

Ki \\\x - a,; II - nl 



X* 



spm spm III 



(36) 



yielding the desired ^i-type cost function Xll^i = (Xll^i -^i)^- sequentially perform the iterations 
(l34l) and (l36l ). starting with A = ^lm> until ||x'^"''^ — x'^H is within some prescribed tolerance e. In our 
simulated scenarios on the order of 3-10 iterations are needed for e = 10~^. This method is denoted by 
SL-li AD. 

2) Non-iterative formulation (SL-ii MD): For a non-iterative solution of (l30l) we start from the 
equivalent formulation (|33] ). with Ai, . . . , Am included as optimization variables through the weighting 
matrix H, and introduce an epigrapth variable ti for each term contributing to tr(-) in the cost function 

minimize tl^, 
U,A,t 

subject to ef (A + RU)'^H(A + RU)ei < U 

\\ui\\ = 1, Xi> 0, l^A = 1, 



(37) 
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where t = ti ... t„ and is the standard coordinate vector with 1 in the i-th position and zeros 
elsewhere. As in [9| we invoke the matrix inversion lemma to express (l32l ) as the limiting case of 
(positive semidefinite) H = lima— >oo(-^ + f'"lmlm)~^' which is more amenable to analytic manipulations 
in optimization problems. In practice we take o" as a sufficiently large constant. Using Schur complements 
the inequality constraint in (|37] ) may be successively written as 

ti ef(A + RU)^ 

(A + RU)ei H-i 

ti{A + almlJn) - (A + RU)eief (A + RU)^ h 0. 



>- 



(38) 
(39) 



The last inequality is bilinear in ti and Ai, . . . , Am, and we linearize it by replacing the optimization 
variable A with a new /3j = tiX. Now, the /3j can be assembled into a matrix 



At, 



(40) 



f3i ... /3, 

which, as shown above, should have rank 1 and satisfy (3ij > 0, l^P = However, the rank-1 constraint 
for (3 cannot be directly imposed in convex formulations, and we resort to a common technique to 
indirectly induce low rank in optimal solutions by adding to the cost function the (scaled) nuclear norm 

Mn. 

Regarding the second term on the left-hand side of (|39l ). we first note that 

1 

Vi 



(A + RU)ei 







1 






Aej R 








Oii R 











(41) 



where Oj and Vi denote the i-th columns of matrices A and U, respectively. Now, consider the following 
variable, obtained from the stacked rotation vectors that make up U, 



W 



1 

vec(U^) 



1 vec(U' 



T\T 



1 



ui uiuf 



W„ 



(42) 



Further, let Jj denote the set of row indices that extracts the elements of [1 vf]'-^ in (|4TI ) from the first 
column of W. Then, the dyad below is readily obtained by selecting the submatrix formed from the li 
rows and columns of W 

1 



Vi 



1 v! 



(43) 
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and this carries over to ( [39l ) through (|4T]) . which can therefore be written in terms of submatrix W/./^. 
The positive semidefinite matrix W will replace U as an optimization variable, retaining the constraints 
along the diagonal blocks in (|42]) . namely, tr(Wjj) = 1. Finally, we obtain the full convex relaxation of 
(l37l) by combining all the above elements and dropping the rank-1 constraint for W that is implied by 
(021) 

minimize tl„ + /x||/3||7v 
W,/3,t 

af] (44) 
R 



subject to diag(/3i) + Ualmlm ^ 



OLi R 



W ^ 0, wii = 1, tr(Wii) = 1, > 0, 1^/3 = t. 
This reference formulation for SL-£i in multiple dimensions is denoted by SL-^i MD. 

3) Simplified non-iterative formulation (SL-ii SD): Our simulation results suggest that in most sce- 
narios the accuracy of the solution obtained from (l44l ) is nearly identical to that of a simplified for- 
mulation where a single epigraph variable, t, is used. Referring to dTT] ). we now minimize tr(tl„) or, 
equivalently, t, and replace the first constraint for all i = l,...,n with the single matrix inequality 
(A + RU)'^H(A + RU) ^ tin. Applying Schur complements as in (|38])-(l39l) yields 



A R 



In 

u 



A^ 
R 



^0, 



(45) 



and again we replace variable A with f3 = t\ such that Pi > 0, l^P = t. Now, however, there is no 
need to assemble a matrix as in ( |40b and to include its nuclear norm as a penalization term in the cost 
function. Finally, to obtain a convex relaxation we replace U with the new variable 

In 



w 



In 
U 



Wii 

U 



(46) 



and drop the rank-n constraint on W that follows from (1461 ). The simplified SDP formulation for SL-^i 
in multiple dimensions, denoted by SL-^i SD, is given by 

minimize t 
W,p,t 

(47) 



subject to diag(/3) + talm^m — 
Who, Wn = In, 



A R 



W 



A^ 
R^ 

1, ft > 0, ll(3 = t. 
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Note that the optimization variables W and /3 in (l44l ) have size {mn + 1) x {mn + 1) and m x n, 
respectively, whereas the corresponding sizes in (l47l ) are only (m + n) x (m + n) and mxl. For ambient 
dimension n = 2 or 3 and for m w 5 anchors used in our simulations problem ( |47l) has considerably 
fewer variables than (l44l) . and the gap increases as m and n grow. 

Given the configuration for variable W in both non-iterative formulations of SL-^i (l42l ). (l46l ). the 
required elements of the rotation vectors that make up U can be obtained from the rightmost (block) 
column of W or by factorizing submatrices along the block diagonal. The former approach is usually 
more accurate 1131 . 



In this section the tightness and the accuracy of our source localization algorithms are tested in 2D and 
3D scenarios, and under various noise assumptions. Results are benchmarked against another relaxation- 
based method proposed in fSl, denoted below as SDR (as in 0), and with the Squared Range LS (SR-LS) 
approach of 141 . While the latter does not resort to relaxation, but rather directly optimizes the source 
coordinates using an iterative root-finding procedure, we take its performance as representative of the 
current state of the art in optimization-based source localization. 

In each reported simulation we performed M Monte Carlo runs, where in each run the source and 
anchor locations were randomly generated from a uniform distribution over a square or cube whose sides 
are [—10, 10]. The observed ranges, corrupted by i.i.d. noise, were generated as described in Section HIl 
under appropriate noise probability densities. The tables list Root Mean-Square Errors (RMSE), computed 



V I? ^i=i ~ a^iP, where Xj and Xj denote the actual and estimated source positions in the i-th 
Monte Carlo run, respectively. For ease of reference, the best result among all algorithms tested for any 
given setup will often be shown in boldface. 

Convexity and tightness of SLCP: In this example we characterize the accuracy of the convex relaxation 
used in SLCP and compare its performance to that of the SDR algorithm of HI. Range measurements 
to a variable number of randomly placed anchors were generated as indicated above over M = 1000 
Monte Carlo runs, and corrupted by white Gaussian noise. 

First, we estimate how often the constraint set S (H), which appears in our formulation of the source 
localization problem prior to relaxation (|7]), is convex along its "upper right" boundary where the optimal 
solution lies. As discussed in Section IIII-A[ when this property holds the relaxed solution $ obtained 
by SLCP ^ will have rank 1 and can be factorized to yield the optimal point for the non-relaxed 
problem dV) on the boundary of S. We empirically assess convexity of S by tracing the boundary of the 



V. Numerical Results 
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(partially hypothesized) convex hull T (fTTb and searching for line segments that delimit regions where the 
boundaries of S and T depart due to local concavity of S. Specifically, we solve the support hyperplane 
problem (fT2l ) for a grid of angles < /? < | and detect the presence of a line segment when the distance 
between the intersection points (u(/3), f (/?)) for two consecutive angles (3 exceeds a threshold. For a 
noise standard deviation dgaussian = 10~^, S passed the convexity test in 80% of runs for three anchors. 
The percentage increased to 84% for five anchors, in line with our reasoning in Section IIII-AI that S is 
more likely to be convex as the number of anchors increases. 

Next, we compare the RMSEs of SDR and SLCP. As in |5| we provide results for all Monte Carlo runs 
(denoted by SDR, SLCP) and also for so-called tight runs (denoted by SDRt, SLCPt) where the solution 
for the relaxed localization problem is close to having rank 1, as desired for subsequent factorization to 
obtain the actual source coordinates. We consider a solution matrix to be tight when the ratio between 
its first and second eigenvalues is at least 10^. Table Hill lists the RMSEs and the number of tight runs 
(-^SDR> -^SLCp) over 1000 trials for five anchors and Gaussian noise standard deviations of 1, 10~^, 10^^, 
and 10^'^. SLCP is clearly superior over the full set of trials, but the gap to SDR closes in the subset of 
tight runs, indicating that the advantage is mostly due to a much higher probability of its solution having 
near rank 1. Even for the highest noise power, where the number of tight runs in both algorithms is 
comparable, the ratio of first to second eiganvalues is usually higher in SLCP, leading to lower RMSE. 

TABLE IIL Source localization accuracy for relaxation-based methods (RMSEs listed for total and tight 
runs). 



^ gaiissian 


A'sDR 


A'sLCP 


SDR |8| 


SDRt |8| 


SLCP 


SLCPt 




490 


921 


0.0045 


0.0014 


0.0020 


0.0015 




444 


815 


0.0162 


0.0107 


0.0112 


0.0108 


10-1 


478 


527 


0.1503 


0.0960 


0.1207 


0.0959 


1 


538 


526 


1.6070 


1.1885 


1.2169 


1.1885 



Under the same simulation setup as above, but using only three anchors, we test the alternative search- 
based method described in Section ITlI-B I to obtain the vector of rotation factors from the relaxed solution 
matrix of SLCP, Improvements in total RMSE are under 1% for all noise variances using 2 x 10^ 
grid points on the interval [0, 2it) to evaluate ([TT] ). Foremost, this suggests that rank-1 factorization by 
SVD, which we adopt as our technique of choice to efficiently extract rotation factors, yields results that 
are indeed very close to the best possible strategy for finding 6. 
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Localization in 2D and 3D under Gaussian noise: In the remaining simulations our algorithms are 
benchmarked against SR-LS H, whose global performance exceeds that of SDR IH because it directly 
optimizes over source locations and does not experience degradations related to tightness of the solutions 
in the same way that SDR does. 

Measurements for the 2D case were generated as in the convexity/tightness assessment, using five 
anchors and M = 200 Monte Carlo runs. Table |IV] Usts the RMSEs for SR-LS, SLCP © and its 
multidimensional counterpart SLNN (l27l) . and also the algorithms for Laplacian noise, namely, the 
complex formulation of SL-^i in [31, and its multidimensional counterpart SL-^i AD ( [34l) . ( [36l ). All 

TABLE IV: 2D source localization under Gaussian noise (5 anchors, 200 Monte Carlo runs). RMSEs are 
given for complex (SLCP, SL-£i) and real (SLNN, SL-^i AD) formulations. 



'^gaussian 


SR-LS |4| 


SLCP 


SL-£i |3 | 


SLNN 


SL-£i AD 


10"^ 


0.0032 


0.0023 


0.0014 


0.0023 


0.0057 




0.0138 


0.0109 


0.0136 


0.0113 


0.0133 


10-1 


0.1406 


0.1037 


0.1118 


0.1097 


0.1249 


1 


1.4947 


1.3249 


1.4536 


1.3580 


1.4593 



algorithms outperform SR-LS, which squares measurements (p = 2, g = 2 in (dJ) and thus becomes 
more sensitive to the presence of (Gaussian) noise in range measurements. The fact that SLCP/SLNN 
achieve the best results is not surprising, as these algorithms actually maximize a Gaussian likelihood 
function. Interestingly, SLCP attains slighly lower errors than SLNN, even though the same cost function 
and similar steps are used in the derivation of both algorithms. However, the impact of relaxing the 
rank constraint in the optimization variable to obtain a semidefinite program is not necessarily the same, 
which could explain the observed differences in performance. Similar comments apply to SL-£i and 
SL-£i AD, although the algorithmic differences between the complex and real formulations are larger 
than for SLCP/SLNN. 

Results for 3D source localization, for which the complex formulations cannot be used, are given in 
Table IV] Again, both SLNN ans SL-^i outperform SR-LS, the former having lower RMSE as its cost 
function is matched to the noise statistics. Source locaUzation with the same number of anchors is a less 
constrained problem in 3D than it is in 2D, resulting in higher RMSEs. Regarding the three variants of 
SL-£i (see Section |IV-Ab . there is no well-defined trend on their relative performance. Note also how 
the RMSE of the non-adaptive simplified formulation SL-^i SD (|47] | is quite close to that of the general 
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TABLE V: 3D source localization under Gaussian noise (5 anchors, 200 Monte Carlo runs). 



^ gaussian 


SR-LS EJ 


SLNN 


SL-^i 








AD 


MD 


SD 


10"^ 


0.0040 


0.0036 


0.0038 


0.0038 


0.0040 


10"^ 


0.0295 


0.0274 


0.0285 


0.0290 


0.0292 


10-^ 


0.2612 


0.2290 


0.2376 


0.2401 


0.2390 


1 


3.3279 


2.7431 


2.9492 


2.8748 


2.8801 



formulation SL-^i MD (l44b . even outperforming it, on average, for one of the noise powers. 

Localization in 2D and 3D in the presence of outliers: The same setup for Gaussian noise is adopted 
here, except that ranges are contaminated either by Laplacian noise, or by what we designate as selective 
Gaussian noise. Range measurements for the latter are created as rj = ||x — aj|| + li^j + |e|, where Wi is 
a Gaussian noise term with dgaussian = 0.04 that is present in all observations and e is also a Gaussian 
disturbance, but with higher standard deviation doutiier G [0.3, 1.5], that contaminates only one measured 
range (i.e., e = for all other observations). This statistical model is less tractable than the Laplacian 
noise model, but we include it in some of our simulations as it more realistically reflects how outliers 
occur in real ranging systems. 

Tables IVTl and I VIll list RMSEs for 2D and 3D source localization under both outlier generation models. 
When compared with Tables |IV] and |V] for the Gaussian case, the most striking difference is that the 
variants of SL-^i, designed for Laplacian noise, now outperform both SR-LS and SLCP/SLNN. As in 
the 2D Gaussian case, the complex formulation SL-^i attains lower errors than the real formulation 
SL-£i AD. In 3D (Table IVIII ) the real variants of SL-^i developed in Section IIV-AI provide significant 
RMSE reductions, on the order of 10%, over those of the Gaussian algorithm SLNN. Interestingly, this 
conclusion still holds for the selective Gaussian case, whose outlier generation model is not matched to 
the Laplacian assumption underlying SL-£i. Similarly to the Gaussian case of Table |Vl the three variants 
of SL-^i exhibit similar performance, here with a slight advantage of the alternating direction algorithm 
SL-£i AD. 

VI. Conclusion 

We have proposed SLNN as an extention to 3 and higher dimensions of the ML-based source local- 
ization approach developed in 13 ] by formulating it as an optimization problem using nuclear norms and 
SDR. Similarly, we also extended to higher dimensions the 2D SL-^i localization algorithm for Laplacian 
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TABLE VI: 2D source localization in the presence of outlier range measurements (5 anchors, 200 Monte 
Carlo runs). 



(a) Laplacian noise 



^ laplacian 


SR-LS (2 


SLCP 


SL-^i 1311 


SLNN 


SL-^i AD 


0.2 


1.1398 


1.0839 


0.9240 


1.0770 


1.0687 


0.4 


1.9031 


1.8585 


1.4546 


1.8908 


1.7996 


0.8 


3.1543 


3.1143 


3.0344 


3.0812 


3.0798 


(b) Selective Gaussian noise 


(^outlier 


SR-LS |4| 


SLCP 


SL-£i 13 1 


SLNN 


SL-^i AD 


0.5 


0.2983 


0.2448 


0.1849 


0.2556 


0.2337 


1.0 


0.4662 


0.4561 


0.2508 


0.4516 


0.3714 


1.5 


1.2419 


1.1640 


1.0542 


1.2389 


1.2157 



noise developed in |9|. Our simulation results show that the proposed algorithms provide very accurate 
results compared to other optimization-based localization methods that operate on range measurements, 
although their performance in 2D is not quite as good as that of the complex formulations developed in 
|[3l . ||9l- In 3D scenarios with Gaussian noise SLNN delivered solutions that were about 5% more accurate 
than those of SL-^i, whereas in the presence of outlier range measurements the situation was reversed 
and SL-£i proved to be about 5-10% more accurate under either Laplacian or selective Gaussian models. 
We developed both iterative (SL-£i AD) and non-iterative (SL-^i MD/SD) 3D extensions of SL-^i, which 
exhibited comparable performance, with a slight advantage of SL-^i AD. Complexity considerations (e.g., 
computational load, maximum admissible problem size) will then play an important role when selecting 
one of those algorithms for a particular application. 

We have carried out an analysis of the geometry of our 2D formulation for ML localization under 
Gaussian noise (SLCP), and found that the high probability that a certain portion of the (outer) border of 
its constraint set is convex justifies the observed strong tightness of our relaxation. Our simulation results 
for random anchor configurations indicate that another well-known SDR relaxation for the same problem 
has a significantly higher chance of yielding optimal solutions that do not have the necessary properties 
(unit rank) to accurately recover source positions. Regarding the extraction of spatial coordinates from 
the positive semidefinite matrix computed by SLCP, we examined a search-based alternative to standard 
rank-1 factorization using the SVD. This strategy is feasible for the practically important case of range- 
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TABLE VII: 3D source localization in the presence of outlier range measurements (5 anchors, 200 Monte 
Carlo runs). 



(a) Laplacian noise 



*^laplacian 


SR-LS (3 


SLNN 


SL-£i 








AD 


MD 


SD 


0.25 


1.3619 


1.3577 


1.2113 


1.2097 


1.1776 


0.5 


2.6514 


2.5719 


2.3236 


2.4704 


2.3651 


0.75 


3.5968 


3.5070 


3.1265 


3.2173 


3.1987 


(b) Selective Gaussian noise 


^outlier 


SR-LS 14J 


SLNN 


SL-^i 








AD 


MD 


SD 


0.3 


0.3930 


0.3237 


0.3033 


0.3035 


0.3033 


0.6 


0.9756 


0.9154 


0.8786 


0.9084 


0.9001 


0.9 


1.5934 


1.5502 


1.2755 


1.2857 


1.3022 



based localization using three anchors, but was found to yield only minor improvements relative to the 
SVD-based factorization. 

Appendix A 
Proof of Lemma [T] 

Our proof of Lemma [1] relies on a result, interesting in its own right, that characterizes the convex 
hull of the set of 3 x 3 rank-1 matrices built from complex vectors with unit-magnitude components. 

Lemma 2. Let 

A = {ee^ : 6> G C^ 1^,1 = 1} , (48) 
S = G C^''^ : * ^ 0, (/>ii = 1} . (49) 

then B = co{A). 

Proof: co(^) C S is straightforward since B is convex and A C B. For the reverse direction 
co{A) D B our goal is to find, for every ^ £ B, matrices @i ^ A and nonzero scalars Aj > 0, with 
Y,i Ai = 1, such that ^ = Y.i 
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Note that both A and B are invariant under the (unitary) similarity operation 

M ^ PMP^, 



(50) 



where P is the product of a permutation and a diagonal unitary matrix. In other words, we can simul- 
taneously permute rows and columns and multiply the i-th row and z-th column by a unit-magnitude 
complex number. Thus, we can assume without loss of generality that * is of the form 

1 a b 

0<a<b, z eC. 



^ = a I z 
b z 1 

Since $ ^ 0, we must have a < 1, & < 1, |z| < 1 and 

-2 b'^ -\z\^ + 2abRe{z}, 



(51) 



< 1*1 = 1 - 



which, for z = X + jy, reads 



{x-abf + y'^ < (1 -a2)(l -62). 



(52) 



(53) 



For fixed a, b this inequality describes a circle (with interior) in the (x, y) plane, centered on {ab, 0). 
Since any point in the interior of a circle can be written as a convex combination of two points on its 
boundary, we can assume that we have equality in (l53l ). Thus, from now on we assume 



z = a 



6+ V(l-o2)(l-62)e 



(54) 



We now complete the proof by expressing such $ as a convex combination of two matrices from A. For 
given < a < 5 < 1 and ip G [0, 2tt[ we want to find a, (5,^,5^ [0, 27r[, and < A < 1 such that 



(55) 





1 


a 


b 




1 






1 




$ = 


a 


1 


z* 


= A 






+ (1-A) 








b 


z 


1 















We thus have 



a = Ae^° + (1 - A)e^'^, b = Xe^^ + (1 - X)e^\ 



From the first two relations we get 

1-A ' 



- Xe^^ 



(56) 
(57) 

(58) 
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and replacing these in the third relation yields, after simple manipulations, 

z = ab + (e"^" - a) {e^^ - b) . (59) 
1 — A 

Before proceeding, we state and prove a useful lemma from elementary geometry: 

Lemma 3. Referring to Figure if A is a point inside a unit circle whose distance to the center is a, 
RS is any line through A, and PQ is a diameter through A, then 

AR- AS = AP ■ AQ = {I- a)(l + a) = 1 - a^. (60) 



Proof Triangles APR and AQS, depicted in Figure QbJ are similar, hence 

AP AR 



AS AQ- 



We use the lemma above with parameters as depicted in Figure |4l From A = XR + (1 — A)^ we have 
4f = and by Lemma |3]Ai? ■AS = l-a'^, hence 



AR=^lj^{l-a^), = a + ^^-j^{l-a^)e^^\ (62) 

Similarly, with A = b, R = e^^, S = e^^, and instead of ipi, we have 



eif^ = b+x/^-j^il-b^)e^^\ (63) 
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R = e^ 




Fig. 4: Application of Lemma |3] to a convex combination on the unit circle. 



Substituting (|63]l back in ^ yields 

z = ab + ^y{l-a'^){l-b'^)e^^^'-'^'\ (64) 



which has the same form as (I541 ). obtained from the positive semidefinite condition for matrix $ in (I5TI ). 

We now argue that letting angle a go from to 27r is equivalent to letting ipi cover an interval of 
length 27r as well (Figure |4ll. Fixing ipi, and consequently a, the two relations in (l56l ). together with an 
arbitrary requirement that Imje-'^} > 0, fix the value 



of (3, 7, 6, A, and, in particular, of ip2. Thus, 



= /(vi) is a continuous function of (^i. 

When ipi = 0, (/?2 has a certain value, say, eq € [0, vr] (it can be computed, but is not needed in this 
proof). For = vr it is straightforward to see that ip2 = tt — eo, and for ipi = 2it it is again eo- In 
particular the continuous function (p2 — <fi takes values from eq — = Eq to eq — 2it, i.e., modulo 2it 
it takes all values in [0, 2-k[. Thus, for any given angle 99 in (l54l ). let ipi be such that f{^i) — (pi = (p, 
modulo 2-K. Then, the corresponding a, /?, 7, 5, and A, as explained above, give the desired decomposition 
(ESJ). ■ 

We now proceed and prove Lemma [T] under the assumption of Lemma [2j thus tacitly assuming m = 3. 
Note that the proof is valid for arbitrary c, r in ([8]) and (ITTI ). i.e., it does not require that the structure 
for these vectors defined in ^ be taken into account. 

^Equivalently, note that fixing fully defines the geometrical construction shown in Figure |4l and thus fixes the values of 7 
and A. Then, A fully defines the corresponding construction for A = h if, in addition, Imje^''} > is specified, and thus fixes 
the values of /3, 5, and tp2- 
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Proof: We rewrite sets 5 in (l8]l and T in (ITTT ) using the notation (1481) 

5 = { (c^0c, r^0r) : G ^} , (65) 
T = { (c^*c, r^*r) : * G co(^)} . (66) 
Obviously S CT. Now let a G [0, ^] and define 



(ui, vi) = arg max ((cos a, sin a), (n, f)). (67) 



We wish to show that 



((cosa, sina), {ui, vi)) = max ((cosa, sina), (u, v)), (68) 

(m,'«)g5 

so that the inner product over S attains the same maximum value as over the larger set T, and the support 
hyperplanes with normal (cosa, sina) thus coincide for the two sets. It is enough to prove that there 
exists {u', v') G S that attains the left-hand side of (|68] ). 
We may write $i G co(^) which maximizes ( [671 ) as 

^i = Y, ^^^i^f^ Ai > 0, J] Ai = 1, \eik\ = 1, (69) 



hence 



((cosa, sina), (ui, vi)) 

= ("v/cos ac)^ cos a c) + (\/sina r)"^$i(\/sina r) 



(70) 



= Y,\i (p^^i^f P + q^^i^f q) 

i 

= ^X, (Ip^e.p + |q^0,,|2) . 

i 

Let io be the index where the last summation attains its maximum value. Then 

((cosa, sina), (m, vi)) < {p^Oif + {ci^Oif = ((cosa, sina), {c^Oi^O^^c, r^0i„6>,f r)), (71) 
which completes the proof because the second argument in the inner product is an element of S. ■ 

Appendix B 
Analysis of SLNN 

Solution of the inner subproblem (|23] ).' For any optimization problem with differentiable objective 
and constraint functions for which strong duality holds, any set of primal and dual optimal points must 
satisfy the KKT conditions |[T6l . We define the Lagrangian of (|23] ) with dual variable A as 

L(V, A) = tr(C^UV) + tr(A^(V^V - I„)) . (72) 
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The first-order KKT conditions are given by 

Vvi(V,A) = U^C + V(A + A^) = 0, (73) 
V^V = I„, (74) 

where (173] ) is obtained by setting to zero the gradienj^ of ( 1721 ) with respect to V, whereas ( 1741 ) is the 
original orthogonality constraint in (|23] ). 

Premultiplying (1731 ) with V^, taking the trace (i.e., taking the inner product with V), and using (1741 ) 
yields the optimal value for the cost function 

tr(C^UV) = tr(V^U^C) = -tr(A + A^). (75) 

But from U^C = — V(A + A^) in (1731 ) we can square both sides to get 

C^UU^C = (A + K^f. (76) 

Hence, among candidate optimal points satisfying the KKT system, the cost function can be made as 
small as possible by choosing A+A^ in (1751 ) as a positive semidefinite matrix square root of the left-hand 
side of (1761) . Replacing this in (1751 ) gives the final optimal cost 

tr(C^UV) = -tr((C^UU^C)2) = -||C^U||iv (77) 

■ 

Interestingly, we point out that the more usual Frobenius norm solves the following relaxed version of 
the inner subproblem (l23l) 

minimize tr(C^UV) = (V, U^C) 

V (78) 
subject to tr(V^V) = ||V|||, = n, 

which is easily verified by writing the KKT system based on the Lagrange function tr(C^UV) + 
A(tr(V^V) - n), 

U'^C + 2AV = 0, tr(V^V) = n, (79) 

whose solution at the minimum is 

r- u^c IIU^CIIf 

^We use the standard results J^tr(A^X) = A and J^tr(XBX^) = X(B + B^) QtI, fT9l . 



November 30, 2011 



DRAFT 



28 



with optimal cost — -v/n||U^C||ir. The minimum cost within the expanded domain of this relaxed 
subproblem will at least be as low as that of (1231) . hence ||U^C||Ar < -/n||U^C||i;'. On the other 
hand, 



IIU^CII^ = ^{^a.f > = IIU^CII^, (81) 

where cjj denotes the i-th singular value of U^C. Combining the two inequalities we have the bounds 

||U'^C||f < IIU'^CIItv < V^llU'^CIIi.. (82) 



Proof of equivalence between ( |26l ) and ( [27l ).- We first rewrite ( |27l ) replacing the linear matrix 
inequality with an equivalent Schur complement 

maximize 2tr(Z) + ^tr(rr-^W) 
W,Z 

(83) 

subject to W ^ 0, Wii = 1 

7? < C'^WC, Z ^ 0. 
Let p\ and P2 be the optimal values of problems (|26l ) and ( [83] ). respectively. 

Choose a feasible point (Z, W) for the second problem, such that ^ Z^ ^ C^WC. This impliej^ 
Z < (C^WC)i, hence the values of the two objective functions satisfy 

2tr(Z) + — tr(rr'^W) < 2tr((C^WC)^) + — tr(rr'^W). (84) 

m ^ m 

In particular, choosing for (Z, W) the unique maximizer of (l83l) . inequality (l84l) asserts that > p^. 

For the converse choose a feasible point W for the first problem and consider the eigendecomposition 
C^WC = QAQ'^. Now set Z = QAiQ'^, so that Z^ = QAQ'^ = C'^WC, and (W, Z) is therefore 
feasible for ( [83l ). For both problems the value of the cost function is 

2 tr( A 2 ) + — tr(rr^ W) . (85) 

m 

In particular, choosing for W the maximizer of (l26l ) the construction for Z yields a feasible point (W, Z) 
for (l83l) where the objective function equals p\. Therefore p\ < p*2, and coupUng this with the converse 
inequality above we conclude that p\ = p\ and the two problems are equivalent. ■ 

"•a >- B ^ A5 >- Bi >- l20l. 
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